www.gusucode.com > Matlab在化学工程中的应用 > Matlab在化学工程中的应用/实用化工计算机模拟－Matlab在化学工程中的应用/Examples/Chapter 3/RTD1.m

    function RTD1
% t1,t2 ----- Time, min
% C1,C2 ----- Concentration, (kmol/m^3)*1e3
% dt    ----- Mean residence time, min
% ss    ----- Variance (sigma square), min^2
% DL2uL1,DL2uL2    ----- Dispersion number
% 
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/05/02 $

clear all
clc

t = [0:20:200];
C = [0, 0, 0, 0, 0.4, 5.5, 16.2, 11.1, 1.7, 0.1, 0];

% 观察样条插值/拟合效果
% By spline():
plot(t,C,'o')
sp = spline(t,C)
hold on
fnplt(sp,'b-')
xlabel('Time (s)')
ylabel('C (kmol/m^3)×10^3')
hold off

 % By spaps():
figure
plot(t,C,'o')
sp = spaps(t,C,1)       
hold on
fnplt(sp,'b-')
xlabel('Time (s)')
ylabel('C (kmol/m^3)×10^3')
hold off

% Method 1: 直接用式（1）、（2）和（3）计算
t0 = 0;
tf = t(end);
IC =  quadl(@Func,t0,tf,[],[],sp); 
ICt = quadl(@Func1,t0,tf,[],[],sp); 
ICt2 = quadl(@Func2,t0,tf,[],[],sp);
tm1 = ICt/IC
ss1 = ICt2/IC - tm1^2
DL2uL1 = ss1/tm1^2/2

% Method 2: 用式（4）、（5）和（3）进行计算
tm2 = sum(C.*t)/sum(C)
ss2 = sum(C.*t.^2)/sum(C) - tm2^2
DL2uL2 = ss2/tm2^2/2

% ------------------------------------------------------------------
function y = Func(x,sp)     % f= C
y = fnval(sp,x);            % 定义被积函数

% ------------------------------------------------------------------
function y = Func1(x,sp)    % f= Ct
y = fnval(sp,x).*x;

% ------------------------------------------------------------------
function y = Func2(x,sp)    % f= Ct^2
y = fnval(sp,x).*x.^2;